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Abstract 

Context. Detailed modeling of the high-energy emission from gamma-ray binaries has been propounded as a path to pulsar wind 
physics. 

Aims. Fulfilling this ambition requires a coherent model of the flow and its emission in the region where the pulsar wind interacts 
with the stellar wind of its companion. 

Methods. We developed a code that follows the evolution and emission of electrons in the shocked pulsar wind based on inputs from 
a relativistic hydrodynamical simulation. The code is used to model the well-documented spectral energy distribution and orbital 
modulations from LS 5039. 

Results. The pulsar wind is fully confined by a bow shock and a back shock. The particles are distributed into a narrow Maxwellian, 
emitting mostly GeV photons, and a power law radiating very efficiently over a broad energy range from X-rays to TeV gamma rays. 
Most of the emission arises from the apex of the bow shock. Doppler boosting shapes the X-ray and VHE lightcurves, constraining 
the system inclination to i ~ 35°. There is a tension between the hard VHE spectrum and the level of X-ray to MeV emission, which 
requires differing magnetic field intensities that are hard to achieve with a constant magnetisation cr and Lorentz factor Y p of the 
pulsar wind. Our best compromise implies cr » 1 and T p » 5 x 10 3 , respectively higher and lower than the typical values in pulsar 
wind nebulae. 

Conclusions. The high value of cr derived here, where the wind is confined close to the pulsar, supports the classical picture that has 
pulsar winds highly magnetised at launch. However, such magnetisations will require further investigations to be based on relativistic 
MHD simulations. 

Key words, radiation mechanisms: non-thermal — stars: individual (LS 5039) — stars: winds, outflows — gamma rays: general — 
X-rays: binaries — methods: numerical 


1. Introduction 


Gamma-ray binaries are composed of a massive star in orbit with 
a compact object and characterized by dominant radiative output 
in the gamma-ray (<; MeV) range (see[Dubus 2013 for a review). 
The compact object is widely thought to be a rotation-powered 
pulsar, although this remains to be proven for most systems. The 
interaction of the pulsar wind with the massive star wind ends up 
dissipating part of the pulsar’s rotation power through particle 
acceleration. Gamma-ray binaries offer an opportunity to study 
the processes involved in pulsar wind nebulae on much smaller 
scales, notably how an e + e~ relativistic wind is launched from 
a rotating magnetosphere and how its energy is released at the 
termination shock. 

Observations of gamma-ray binaries show flux variations 
tied to the orbital phase in most of them. Theoretical models 
have focused on relating these variations to changes in the line- 
of-sight and/or the conditions at the termination shock as the 
pulsar follows its eccentric orbit. Amongst the different systems, 
LS 5039 constitutes a useful testbed because of its regular, well- 
documented modulations in X-rays (Takah ashi et al.|2 009 ), low 
energy (LE, 1-100 MeV,|Collmar & Zhang|2014| l, high energy 
(HE, 0.1 - 100 GeV, Abdo et al. 2009) and very high energy 


(VHE , > 100 GeV, |A taronian e t al.|20Q6[ ) gamma rays. The X 


ray, LE and VHE modulations are in phase, with a peak at infe¬ 
rior conjunction (when the compact object passes in front of the 
massive star as seen by the observer) and a minimum close to su¬ 
perior conjunction (compact object behind the star). In contrast, 
the HE gamma-ray modulation is in anti-phase, peaking at supe¬ 
rior conjunction. LS 5039’s short orbital period of 3.9 days has 
made it easier to establish these modulations than in the other 
systems. The modulations appear stable over time, with no re¬ 
ported change in the orbital lightcurves. 

Nearly all models invoke synchrotron and inverse Compton 
emission from pairs with energies up to several TeV to explain 
the high-energy emission from LS 5039. These processes are 
efficient in the sense that the energy loss timescale is usually 
short compared to the flow timescale for the most energetic pairs: 
a large fraction of the available power may thus end up as high- 
energy radiation (e.g. [Bosch-Ramon & Khangulyan|2009] l. 

The seed photons for inverse Compton emission are provided 
by the massive star. The pairs upscatter these stellar UV pho¬ 
tons to HE and VHE gamma rays. The gamma-ray emission is 
maximum when pairs backscatter the stellar photons in the di¬ 
rection of the observer i.e. at superior conjunction. However, the 
gamma-ray photons must propagate through the binary system 
before reaching the observer. VHE photons are above the thresh¬ 
old for pair production with stellar photons (<; 30 GeV when in- 


1 

















Dubus, Lamberts, Fromang: Modelling gamma-ray binary emission based on relativistic hydrodynamics 


teracting with the kT+ ~ 3 eV blackbody photons from the star), 
hence a large part of the VHE flux can be lost to creating e + e~ 
pairs before the radiation escapes. The yy opacity in LS 5039 is 
minimum close to inferior conjunction and maximum close to 
superior conjunction, explaining why HE and VHE gamma rays 
are anti-correlated although both arise from inverse Compton 
emission of the same seed photons (see §4 in Dubus|[20l3 and 
references therein). 

There are two caveats to this interpretation of the HE and 
VHE gamma-ray modulation. First, a straightforward applica¬ 
tion predicts that none of the VHE flux emitted in the vicinity of 
the compact object should make it through the system at supe¬ 
rior conjunction, whereas observations still detect a faint source. 
One explanation is that there is a contribution to the emission 
from the electromagnetic cascade that occurs as newly-created 
e 1 e pairs emit VHE photons which in turn create pairs etc (e.g. 
|Bednarek|2006[ [Bosch-Ramon et al.|2008[|Cerutti et al.|2010j i. 

Another, non-exclusive, explanation is that the VHE emission 
arises from a larger or more distant region, diluting the effets 
of the yy opacity (e.g. jZabalza et al.|2013) . The second caveat 
is that the HE and VHE spectra are very distinct, with the HE 
spectrum cutting off exponentially at a few GeV ( |Hadasch et al.| 
2012). Clearly, different populations of electrons must be in- 


R, 

d 


4/2 


1 + rj 


1 / 2 ’ 


( 1 ) 


where rj — E/M w v w c, with E the pulsar spindown power, M u 
and v w the stellar wind mass loss rate and velocity. The termi¬ 
nation shock distance R s doubles from periastron to apastron. 
decreasing B oc 1/R S proportionally (e.g. Dubus 2006b; Takata 


& Taam 2009). The changing separation can also affect adiabatic 


cooling of the particles as they are advected away from their ac¬ 
celeration site ( |Takahashi et al.|2009| . However, it is difficult to 
tie such changes to the observed X-ray modulation; the peaks 
and dips at conjunctions suggest that it is due to a geometrical 
effect related to the observer line-of-sight rather than due to in¬ 
trinsic changes in the conditions experienced by the particles. 
One possibility related to the line-of-sight is Doppler boosting. 
If rj is small, the termination shock has the appearance of a bow 
shock facing away from the massive star. At inferior conjunc¬ 
tion the bow shock flows in the general direction of the observer, 
whereas it is flows in the opposite direction at superior conjunc¬ 
tion. If the shocked pulsar wind retains a moderately-relativistic 
bulk motion, the synchrotron emission will be boosted at inferior 
conjunction and deboosted at superior conjunction. Dubus et al. 
( 2010j ) have shown that the effect is significant enough to be able 
to explain the X-ray modulation. 

A more accurate assessment of this scenario requires numeri¬ 
cal simulations. The geometry of the termination shock, the bulk 


velocity of the shocked fluid at each location, the adiabatic losses 
can all be derived from a relativistic hydrodynamic al simulation 
instead of being parametrized as done in previous works. The 
impact of yy absorption, Doppler boosting and particle cool¬ 
ing on the orbital lightcurve can then be quantified properly, 
tightening the constrains on the underlying particles and pul¬ 
sar wind physics. Simulations of gamma-ray binaries have been 
performed by several groups, focusing on the geometry of the 
interaction region. Relativistic hydrodynamical simulations in¬ 
dicate that the material re-accelerates to very high Lorentz fac¬ 
tors in the tail of the bow shock (Bogovalov et al. 2008! and 


that the non-zero thermal pressure in the stellar wind leads to 
smaller opening angles for the bow shock than usually assumed 
( Lamberts et al.|2013) l. Pulsar wind nebulae are thought to have 
a low magnetisation <x at the termination shock (Gaensler & 
Slane 2006), in which case the magnetic field has a negligi¬ 
ble impact on the shocked flow, and indeed Bogovalov et al. 
( |2012] > found little difference, on scales of order of the orbital 
separation, between their relativistic magneto-hydrodynamical 
(RMHD) and relativistic hydrodynamical (RHD) simulations of 
colliding winds in gamma-ray binaries. Other simulation work 
includes [Bosch-Ramon et ak ( 2012) 1, who studied mixing due to 
orbital motion using large scale 2D relativistic hydrodynamical 


volved in the HE and VHE domains. Their origin is uncertain. 
Possible sites that have been considered include: the pulsar mag¬ 
netosphere, the pulsar wind, various locations along the pulsar 
wind termination shock, the stellar wind termination shock. 

The interpretation of the X-ray and LE gamma-ray modula¬ 
tion also requires additional ingredients. Synchrotron emission 
dominates in this range. Its luminosity depends on the density of 
electrons and the magnetic field B, both of which can be affected 
by the changing distance of the termination shock to the pulsar. 
Indeed, with an eccentricity e = 0.35 (Casares e t al.|2003| , the 
orbital separation d in LS 5039 varies from 0.1 AU at periastron 
to 0.2 AU at apastron. If the winds are isotropic and have a con¬ 
stant velocity, the distance R s to the termination shock is set by 
(e.g. [Lebedev & Myasnikov|1990]) 


simulations, Paredes-Fortuny et al. (2015 1 who looked at the im¬ 
pact of a clumpy stellar wind on the shock structure, and Takata 
et al. \2012) who presented two 3D SPH non-relativistic sim¬ 
ulations of the interaction of a pulsar wind with a Be disc and 
wind (covering very large scales with limited spatial resolution). 
Takata et al. ( 2012j > computed some lightcurves and spectra from 
their simulations, but do not take into account particle cooling 
and relativistic effects. A comprehensive approach linking high- 
resolution simulations of the shock region with particle cooling 
and emission (including relativity and anisotropic effects) has 
not been attempted yet. 

Here, we use a relativistic hydrodynamical simulation as the 
basis for such a comprehensive model of the emission from par¬ 
ticles in the shocked pulsar wind, with the aim of explaining 
the spectral energy distribution and the flux orbital modulations 
observed from LS 5039. The simulation (§2) provides the spa¬ 
tial evolution of the density, velocity and internal energy nec¬ 
essary to compute the extension of the emission, the impact of 
Doppler boosting, the importance of adiabatic cooling. This ra¬ 
diative post-processing step is described in §3. The general re¬ 
sults are presented in §4, the more detailed application to LS 
5039 in §5. 


2. Relativistic hydrodynamics 

We perform a 3D simulation of LS 5039 using the RAMSES hy¬ 
drodynamical code ( jTeyssier 2002) . The complete description 
of the extension to special relativity is in Lamberts et al. ( 2013) 1. 
Here, we recall the major aspects of the method. RAMSES uses 
a Cartesian grid and allows Adaptive Mesh Refinement (AMR) 
so as to locally increase the resolution at a reasonable compu¬ 
tational cost. It solves the 3D-RHD equations using an upwind 
second order Godunov method. These equations can be written, 
for an ideal fluid, as a system of conservation equations in the 
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laboratory frame ( jLandau & Lifshitz|l959| 

dD d(Dv k ) 
dt 


d_Mj 

dt 


+ 


dx k 

d(M jVk + p6 j - k ) 
dx k 

dE d(E + p)v k 
dt dx k 


= 0 

(2) 

= 0 

(3) 

= 0, 

(4) 


where the vector of conservative variables is given by 



D j 


' Yp 

U = 

Mj 

= 

Y 2 phvj 


UJ 


X 2 pc 2 h - p, 


(5) 


D is the density, M the momentum density and E the energy 
density in the frame of the laboratory, c is the speed of light, 
the subscripts j, k stand for the dimensions, d jk is the Kronecker 
symbol, h is the specific enthalpy, p is the proper mass density, Vj 
is the fluid three-velocity, p is the gas pressure. The gravitational 
force is neglected since we do not treat the acceleration of the 
winds and the outflow speeds are well above the escape velocity 
of the system. The Lorentz factor is given by 


r = 


y/l - V 2 /c 2 


(6) 


The set of Eqs. [2]4 is completed by an equation of state that de¬ 
scribes the thermodynamics of the fluid. In RHD, the ideal gas 
approximation is inconsistent with the kinetic theory of relativis¬ 
tic gases ( Taub|1948] >. The exact equation of state for relativistic 
fluids involves modified Bessel functions, which lead to addi¬ 
tional computational costs. Therefore we use the approximation 
developed by Mignone et al. ( (2005] ), which differs from the ex¬ 
act solution by less than 4 %. In this case, the specific enthalpy 
is given by 


(7) 



+ 1. 


With the specific internal energy defined as e = /t / | (y - l)pc 2 ], 
this gives a variable equivalent adiabatic index of 


7 = 


h- 1 


pc 2 


( 8 ) 


In our simulation, we use a HLL solver with a second or¬ 
der reconstruction based on a minmod slope limiter. This solver 
limits the development of Kelvin-Helmholtz instabilities in the 
simulation (Lamberts et ak| |2013j l, a deliberate choice that we 
justify in jj |4.1| The size of our cubic simulation box / hox is six 
times the binary separation d. With z the binary axis, the box 
extends along the binary axis from z = -1.5 to z = 4.5 in a 
coordinate system scaled by the orbital separation d. The star is 
at the origin and the pulsar at x = 0,y = 0,z = 1. The jc-axis 
extends from -3 to 3 and the v-axis extends from y = -1.5 to 
y = 4.5, so the binary is not located at an edge of the simulation 
box. Our coarse grid is made of 128 cells and we use four levels 
of refinement. This gives an equivalent resolution of 2048 3 cells. 
The high resolution zone, where the four levels of refinement are 
effectively used, covers a slab of width Ay = Az = 2 centered on 
the pulsar and Az = 1 centered on x — 0. Refinement is based on 
density and Lorentz factor gradients. 

We initialize the winds in a spherical region with a contin¬ 
uous outflow. The spherical regions both have a radius of 0.15 


(in units of d). This is large enough to yield spherical symme¬ 
try, but small enough to avoid an impact on the formation of the 
shock. The winds are updated in the spherical region at every 
timestep. We assume the winds have reached terminal velocity 
at launch. A constant low-density medium initially fills the sim¬ 
ulation box. This medium is cleared out as the winds propagate 
out and interact. After a transition phase corresponding to about 
lOfdyn (where we define fd yn as the time for the pulsar wind to 
reach the back shock located at x * 3), the simulation converges 
to a laminar, stationary state where the position of the shocks 
does not evolve with time anymore i.e. time drops out of Eq. 2- 
4. Our radiative calculation is based on a snapshot in the (y, z) 
plane of the hydrodynamics in this state. The geometry of the 
interaction depends only on the ratio of wind momentum fluxes 
ri (§1). For LS 5039, tj is thought to be - 0.1, implying a reason¬ 
able value for the pulsar spindown power E « 4 x 10 W| erg s 1 


2000 km s (Szostek & 


assuming M w & 10 M Q yr , v H 
|Dubus|201 l[|Zabalza et al.|2011[ >. 

In the simulation, the stellar wind has a mass loss rate of 
M w — 2 x 10 8 M e yr 1 and a velocity of 2000 km s '. It has a 
Mach number Al = 30 at the location of the pulsar. Conventional 
models of pulsar winds invoke ultra-relativistic Lorentz factors 
up to T p = 10 6 . Such high values are well beyond the range that 
can be achieved by standard fluid methods. We set the velocity of 
the pulsar wind to v p = 0.99c, or T p = 7.08. This is high enough 
to capture the relativistic effects in the shocked flow, especially 
near the apex where the shock is perpendicular and the post¬ 
shock flow speed tends to (y - l)c when Y p » 1. In the wings, 
the shocked flow gradually accelerates to a velocity close to the 
initial pre-shock velocity. The spatial scale on which this occurs 
increases with the initial Lorentz factor of the wind ( [Lamberts] 
et al.|2013) . The pulsar mass loss rate is scaled to have p = 0.1, 


so that 


M p = 7] 


M w v w 

TpVp 


(9) 


The flow dynamics and geometry will be correct even if the den¬ 
sity must be scaled. For reference, the pulsar wind power corre¬ 
sponding to our choice of M„. = 2 X 10 8 M 0 yr 1 and // = 0.1 
in the simulation is E = Y p M p c 2 » 7.6 x 10 35 erg s . However, 
note that we can scale the spindown power up or down without 
changing the hydrodynamical simulation as long as M w is scaled 
in proportion to keep rj — 0.1. The classical Mach number of the 
pulsar wind is set to 20, which gives a relativistic Mach number 
Al r ei = MY p - 140. Orbital motion turns the shocked struc¬ 
ture into a spiral of stepsize of order S ^ v s P or b — 4 AU =* 20 
d ( |Lamberts et al.||2012[ l. Our simulation covers a smaller size 
(- 5 d) and we make the assumption that we can neglect or¬ 
bital motion on this scal^j Even if previous simulations indicate 
that this is a reasonable approximation when the weaker wind is 
also the fastest (see e.g. Fig. 6 in [Lamberts et al.||2012 or Fig. 
1 in |Bosch-Ramon et aI7||2012j >, a more realistic model should 
include the orbital motion. We discuss in §6.2 how this might 
change our results. 

The ratio p does not change along the orbit since the winds 
are isotropic and are at their terminal velocity. The dynamical 
timescale corresponds to fdyn ~ 300 s for the typical orbital sep¬ 
aration of 0.2 AU appropriate to LS 5039 (see above). The time 
taken to form the interaction region, on the scales that we con¬ 
sider, is very short compared to the orbital timescale of 3.9 days 


1 Neglecting orbital motion makes the problem 2D axisymmetric 
around the binary axis. RAMSES currently does not allow 2D axisym¬ 
metric calculations so our simulation needed to be 3D to allow a quan¬ 
titative comparison with the observations. 
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Figure 1 . Model geometry showing our choices of axis and an¬ 
gles. The interaction has cylindrical symmetry around the z axis. 
The two thick solid lines in the (x, z) plane represent the pul¬ 
sar wind termination shock and the contact discontinuity with 
the shocked stellar wind. Dashed lines represent the same in the 
(y, z) plane. Particles injected at the termination shock follow 
streamlines in the region in between these two surfaces. 


for LS 5039. The structure has ample time to reach a stationary 
state at each orbital phase. Hence, since the stationary structure 
only depends on //, which is constant with orbital phase, one sim¬ 
ulation is enough for all orbital phases even if the orbital sepa¬ 
ration changes. The velocities, M and wind temperature (Mach 
number) remain constant. However, for each orbital phase we 
rescale the distances to the actual orbital separation i.e. by a fac¬ 
tor d. The density and pressure in a given pixel are thus rescaled 
by a factor l/d 2 . 

By doing a RHD simulation instead of a RMHD simulation, 
we implicitely assume that the magnetic field has no dynamical 
role in the interaction on the spatial scales probed by the simu¬ 
lation. This assumption is supported by the study of Bogovalov 
et al. (|2012|>, who carried out both relativistic hydrodynamical 


and relativistic MHD simulations of interacting winds, on spa¬ 
tial scales comparable to those studied here, and found little dif¬ 
ference between the two when the pulsar magnetisation cr < 0.1. 
Conventional models of pulsar wind nebulae assume cr <s 1 at 
the termination shock to ensure the efficient conversion of the 


wind kinetic energy to particle energy at the shock (Kennel & 


|Coroniti|1984[[Gaensler & Slane|2006||Bucciantini et al.|201 1 ( i 

We come back to this issue in §6.2. 


3. High-energy radiation 

The interaction of stellar winds classically leads to a double 
shock structure separated by a contact discontinuity (Fig. [TJ. 
Non-thermal emission may be expected from the shock asso¬ 
ciated with both winds ( Bednarek|2011) . There is plenty of evi¬ 
dence for efficient high-energy non-thermal emission from pul¬ 
sar wind termination shocks, but rather less from stellar wind ter¬ 
mination shocks with only Eta Carinae detected so far in gamma 
rays (Werner et al .|2013) l. Here, we have assumed that the non- 
thermal emission results exclusively from particles accelerated 
at the pulsar wind termination shock. We have not taken into 
account thermal bremsstrahlung from the shocked stellar wind, 
since this would be best treated by dedicated simulations a la 
Stevens et al. ( 1992[ >, and because the X-ray emission from LS 
5039 appears entirely non-thermal with, as yet, no evidence for 
thermal emission (|Zabalza et al.|20lT|l. 


Particles are randomized and/or accelerated to very high en¬ 
ergies at the termination shock of the pulsar. We wish to fol¬ 
low the evolution of these particles once they are injected in the 
shocked pulsar wind flow. Once their energy distribution at each 
location is known, we can compute the radiation seen at a given 
line-of-sight to the binary. 

We adopted some basic assumptions to make the problem 
tractable. First, the evolution of the particles is decoupled from 
the hydrodynamic al simulation and treated in the test particle 
limit as a post-processing step i.e. radiative cooling does not im¬ 
pact the dynamics of the flow. We come back to this in j |4.2 and 
§5.1. Second, the high-energy particles are assumed to follow 
the flow. Spatial diffusion is expected to remain negligible if the 
Larmor radius remains small compared to the flow spatial scales, 
which is the case here. Third, plasma processes are assumed to 
keep the particle momentum distribution isotropic. Last, the flow 
is assumed to be stationary, simplifying the problem to that of 
following the evolution of particles along streamlines. 

We start by explaining how we calculate the total emission 
based on the particle evolution along streamlines (j p.l| >. We 
then describe how we choose the initial distribution of the in¬ 
jected particles (§3.2), how we compute the particle evolution 
and streamline emission ({j |3.3| >. how we estimate the magnetic 
field (§ |3.4| l, and how we deal with the changing aspect due to 
orbital motion ( §3.5| ). 

3 .1. Total flux 

The total flux from the emission region, measured in the labora¬ 
tory frame at a frequency v, is given by 

F'(v) = ^ J' D 2 bs j(v/D obs )exp(-T v ) dV, (10) 

where the integral is over the volume V of the region (measured 
in the laboratory frame), D is the distance to the source, j is the 
local particle emissivity in the co-moving frame (j is the sum 
of the synchrotron j sync and inverse Compton / lt emissivities), 
and t v is the opacity at frequency v due to pair production as 
the photons emitted in the volume dV travel along the line-of- 
sight to the observer (see s D 0 bs is the Doppler boost to the 
observer 

£>o b s = [r(l-v.e obs )r', (11) 

with v the flow velocity vector and e obs the unit vector in the 
direction of the observer ( §3.5[ ). 

The simulated flow is laminar. To compute this integral, we 
extract from a snapshot N streamlines in the emission region. 
Each streamline starts at the position of the shock and ends 
where it leaves the computational domain. These streamlines 
subdivide the emission region into streamtube^] The integral 
over the volume can be written as an integral over time of the 
particles flowing along each streamtube: 

1 N C u 

L(v)=^^>]J Dl bs j (v/D obs ) exp (~T v )\.Sdt, (12) 

where S is a cross section of the streamtube i and t = tj is the 
time taken by particles flowing along streamline i to leave the 
computational domain, with t = 0 at their injection at the shock. 
Since the flow is stationary, the particle flux along each stream- 
tube 

Ni = Lkv.SI, (13) 

2 Rather than streamtubes, these are actually hollow cones since the 
model has axial symmetry around the binary axis, see Fig. |Tj 
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is constant (n = p/m e in the pulsar wind composed of e + e 
pairs). The total flux from the shocked region becomes 

1 N 

F(v) = — Y.fivW,, (14) 

i= 1 

where we have defined the fluence per particle f of the stream¬ 
line i as 


fi(y) = 



D 2 obs j (v/Dobs) exp (-r y ) 

--- at. 


(15) 


with Ri the gyroradius. We expect £, > 1 for diffusive shock 
acceleration, with £, < 10 corresponding to “extreme accel¬ 
eration” (Khangulyan et al. |2008| l ; £ < 1 may be possible 
for acceleration at magnetic reconnection sites, as proposed to 


explain the gamma-ray flares from the Crab nebula (Cerutti 
et al.|2012]l. Synchrotron losses dominate over inverse Compton 


losses at very high energies (Dubus 2006b). The synchrotron loss 
timescale is 


Tsync — y 



ym e c 2 

| o- T c(J3y) 2 u b 


* 77 



s, 


(18) 


This is an integral over the emission per particle as they flow 
along their streamline. The integral can also be recast as an in¬ 
tegral over position along the streamline since the curvilinear 
distance l is related to time by dl = vdt. We obtain N, numeri¬ 
cally by measuring the particle flux across the shock surface for 
each streamtube in the RHD simulation (Eq. [L3j). 

We make the implicit assumption in Eq. 1 14| that the emis- 
sivity j and the opacity r y do not vary substantially across 
the streamtube for a given position. In practice, taking a suffi¬ 
cient number of streamlines ensures this is achieved. We used 
N — 920 streamlines, evenly spread along the shock surface, 
and we checked that this is more than enough for numerical con¬ 


vergence of Eq. 14 


Alternatively, one could divide up the particle distribution 
into energy bins, recast the evolution equation of the particles 
in each energy bin in conservative form and add them to the set 
solved by RAMSES (e.g. |Reitberger et al. 2014[ >. The advan¬ 
tage is that this can deal with non-stationary flows and mixing. 
However, a major drawback is that in our case the particle cool¬ 
ing time at very high energies can be very short: the particles 
cool on small spatial scales, requiring a very high spatial resolu¬ 
tion (see jj |4.2| below). Another difficulty is that radiative cooling 
does not scale with the orbital separation d as adiabatic cooling 
(d~ 2 and d 1 respectively, see Eq. 23 below). Each orbital phase 
would thus require a full simulation to account for this differ¬ 
ence in cooling spatial length. Each new choice of parameters 
for radiative cooling would also require a new simulation run. 
Treating the particle evolution in a post-processing step makes 
the problem much more tractable, allowing for a wider explo¬ 
ration of the parameter space involved in the flow emission. 

We explain our choices for particle injection at the shock 
first before turning to the computation of the streamline fluence 
(Eq.[T5}inQ 


3.2. Particle injection at the shock 

We distribute the particles as a power-law function of their 
Lorentz factor y: 


dn 

dy t =o 


= Ky~ s 


with y min < y < y max . 


(16) 


where b and t' are the magnetic field and time in the comoving 
frame. Setting r acc < r sync gives 


= / 3e \ m 

\£cr T b) 


5x10' r‘ /2 I — 

b 


1G 


1/2 


(19) 


With y max and K known, we derive y m ; n by writing that the 
energy in the particle distribution is a fraction ij p of the down¬ 
stream specific energy e : 



The numerator is the total energy in the distribution and the de¬ 
nominator is the particle density. This equation implicitly sets 
ymin, which can be found numerically with a Newton-Raphson 
scheme. If y max » y m j n and .v > 2 then the integrals can be 
simplified and solved for y m i n , suc h that 


'■fa 


/,-2) 

1 - ^ 1 

(s-2\ 

1 P 

U-iJ 

1 (r - i)' 


nm e c 2 


( 21 ) 


£p enables us to scale the emission to realistic values since the 
(lab-frame) specific enthalpy of the cold pulsar wind at the shock 
is Y p « 7 in the simulation, which is too low to reproduce the 
very high values of e (equivalently T ; ,) required to fit the obser¬ 
vations. Raising tj p implies that the effective mass loss rate from 
the pulsar wind is decreased by because E/c - p Y p M p c is 
fixed by // for given M w , v w , T p (in other words, the same total 
available energy is distributed amongst fewer particles). 

We have also experimented with a relativistic Maxwellian 
distribution because Fermi-LAT observations of gamma-ray bi¬ 
naries in HE gamma rays require an additional population of 
particles with a narrow range in energy (§1). Moreover, particle- 
in-cell simulations of relativistically-shocked pair plasmas typi¬ 
cally show a prominent Maxwellian distribution of shock-heated 
particles together with the power-law distribution of acceler¬ 
ated particles ( Sironi & Spitkovsky|2009 2011] >. The relativistic 
Maxwellian distribution is 


The particle density at the shock sets the normalisation K of the 
distribution. The two other parameters are y nlm and y max . 

The maximum Lorentz factor y max is set by the balance be¬ 
tween acceleration and radiative losses, since the gyroradius will 
be much smaller than the characteristic size of the acceleration 
region in our case ( |Dubus|2006b| ). We assume that the acceler¬ 
ation timescale in the comoving frame r acc is some multiple of 
the Bohm limit: 

r acc = 2ntjR L /c, (17) 


dn 

dy 


t=o 


Ky 1 exp (-y/y,). 


( 22 ) 


Again, K is derived by imposing that the integral of the distri¬ 
bution scales with the particle density at the shock. The mean 
Lorentz factor of the distribution y, is derived from Eq.[20] which 
gives y, « e nt /3 when y min « y, <K y max . 

We assumed that 2j p and i : do not vary along the shock for 
lack of strong justifications for more sophisticated assumptions. 
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3.3. Streamline emission 

For a stationary flow, the evolution of a particle injected at a 
given location and followed along the associated streamline in 
the comoving frame, is uniquely set by the evolution equation 
(e.g. Del Zanna et al. 2006! ; Mimic a et al.j2009][Porth et al. 2014 
for applications to pulsar wind nebulae or AGN jets) 


1 dy dine 1 
y dt' dt' T syn 


1 


(23) 


where the terms on the right hand side represent adiabatic, syn¬ 
chrotron and inverse Compton losses - the only relevant cooling 
processes here. The elapsed time in the laboratory frame dt is 
related to the proper time in the comoving frame by dt - I dt'. 
The bulk Lorentz factor T and the elapsed time dt are derived 
from the position and velocity along the streamline as calculated 
with RAMSES. The adiabatic loss term is also directly derived 
from the simulation. 

The evolution of the magnetic field along the streamline must 
be known to compute the synchrotron losses r sync (Eq. 18 1 and 


/ mc \ 3 1 

(R*) 

2 
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l h ) ' 

UJ 

exp | 

( E+rtigC 2 \ 

(iwJ 
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£)* is the Doppler boost required to transform the stellar radia¬ 
tion field into the comoving frame. 


a>* = [T(l-v.e*)]- 


(25) 


with e* the unit vector giving the direction from the star to the 
flow element containing the pairs. 

The evolution of the particle distribution can be calcu¬ 
lated semi-analytically when inverse Compton losses are in the 
Thomson regime ( |Begelman & Li|[l992) >. A numerical solution 
is required in our case since stellar photons are upscattered to 
gamma-ray energies in the Klein-Nishina regime (Dubus et al. 
2008). Following |Bosnjak et al.| ( f2009] >, the particle distribution 
of each streamline is discretized in “Lagrangian” bins 


rik 


r with V 

J n d y k 


r>k 


(26) 


The relative number of particles n k /n in each energy bin is con¬ 
served but the bin boundaries [jk, yk+\\ vary along the streamline 
according to the energy loss equation (Eq.|23j). We use 400 en¬ 
ergy bins, initially logarithmically distributed between y min and 
Tmax- For each y/ ( , Eq. [23] is integrated in time steps represent¬ 
ing a small fraction (5%) of the minimum energy loss timescale. 
To ease the computational burden, we stop following the energy 
losses once y becomes lower than 10 3 : we verified that these par¬ 
ticles do not contribute emission in the frequency range we are 
interested in. 


We compute the streamline fluence fi(v) (Eq. 15 i once the 
evolution of the particle distribution along the streamline is 


the associated emissivity j sync . Our choices for the magnetic field 
are explained below in §3.4. The synchrotron emissivity j sync 
is computed using the usual formula involving Bessel functions 
( Rybicki & Lightman|1979| . 

For the inverse Compton losses t; c , we take the massive star 
as the only source of seed photons and compute the electron en¬ 
ergy losses r; c using the Jones ( 1968j > scattering kernel. The star 
is modeled as a blackbody of temperature T* = 39 000K and 
radius R* = 9.3/0. The density «* of photons with an energy 
e* (in units of m e c 2 ) seen in the comoving frame by pairs at a 
distance ii* from the star is, in photons per cm 3 per unit energy. 


(24) 


known. The fluence depends on the location in the flow but also 
on the line-of-sight to the observer because of the relativistic 
boost associated with the bulk motion of the flow D 0 \, s (Eq. 111. 
Unlike the particle evolution calculation, which is symmetric 
around the binary axis and requires only a 2D integration (see 
§ |3.5| l. the emission calculation requires a full 3D integration 
since e 0 b s also varies with the azimuth 0 around the binary axis 

(Fig.fT). 

Synchrotron radiation is isotropic in the comoving frame for 
a random orientation of the magnetic field and an isotropic dis¬ 
tribution of particles, both reasonable assumptions at this stage. 
However, the inverse Compton emission is clearly not isotropic 
in the comoving frame, even for an isotropic distribution, be¬ 
cause the seed photons from the star are anisotropic |^| We fol¬ 
low |Dubus et al.| ( 12010) 1 to take this effect into account when 
computing the upscattered emissivity y) c towards the observer 
line-of-sight. Finally, we also calculate the yy absorption of the 
VHE flux due to pair production with the stellar photons as in 
Dubus l2006al). The yy opacity r rr depends on the path of the 
VHE photon from its emission location to the observer. We ne¬ 
glected the finite size of the star to ease the computational re¬ 
quirements. This approximation appears reasonable here since 
the extension of the VHE emission region would probably re¬ 
duce and smoothe out any effect of the finite size on both the 
inverse Compton emission and the yy opacity. Knowing r y , yj c 


and _/ sync along each streamline enables us to evaluate Eq. 15 


3.4. Magnetic field 

Following the standard description of pulsar winds, the labora¬ 
tory frame magnetic field B p at a distance d p from the pulsar in 
the unshocked wind is given by 


B p tl +cr 
4n 


cr I 4nd 2 c ’ 


where cr is the wind magnetisation 


Bi 


4jtYl 


p n P m e Vp 


(27) 


(28) 


and the subscript p identifies a value in the unshocked wind. 
Pulsar winds are thought to have a low cr so we neglect the in¬ 
fluence of the magnetic field on the flow dynamics (§2). 

The magnetic field is purely toroidal at large distances from 
the magnetosphere of the pulsar, so the shock is perpendicular 
and the field is amplified at the shock by the compression ra¬ 
tio^ = ItlIt p — N/N p , where N = Fn is the laboratory frame 
number density. The compression ratio depends on the adia¬ 
batic index y, which can change along the shock (see Eq. |7J. 
Even if y is constant, the compression ratio changes in a rela¬ 
tivistic shock when the transverse speed is non-negligible. For 
a ultra-high wind Lorentz factor, the shock behaves like a nor¬ 
mal shock with y « 1 /(y - 1) except very far out in the wings, 
where the speed normal to the shock becomes non-relativistic 
and^ = (y + l)/(y- 1). We have taken x = l/Cy- 1) everywhere 
along the shock, so that 


B = 


1 


f - 1 


d 2 p c 


(if?) 


1/2 


(29) 


3 By using the Jones kernel we have assumed that particles (contin¬ 
uously) lose energy isotropically on average. This is consistent with 
anisotropic emission if plasma processes maintain the particle distribu¬ 
tion isotropic. 
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In practice we take B — Ti,/d p , where Ch is a constant, since E 
and cr are not fixed but will be determined by the comparison to 
data. 

Alternatively, the magnetic field could be amplified by 
plasma instabilities to a value that represents a fraction of the 
available internal energy 


b 2 _ , P 
87t y-1’ 


(30) 


where b = B/Y is the magnetic field in the comoving frame, 
again assumed to be toroidal, and (again) is a constant. From 
the conservation of Yli across the shock, the pressure is related 
to the upstream conditions by 


P = 


y - i / r , 


r 


“F T7 1 1 _ tH n P m <- c 


when the upstream pressure is negligibl^] This equation shows 
that p oc {Yd p y 2 , hence, that B oc 1 /d p . Thus, there is little 
difference in practice between the magnetic field B as given by 


Eq. 29 and Eq. 30 the two being exactly proportional when the 
pulsar wind Lorentz factor is high (Y p » 1). 

Since we assume no influence of the magnetic field on the 
flow, the evolution of B beyond the shock is set solely by the in¬ 
duction equation. For a stationary flow, with cylindrical symme¬ 
try around the binary axis, the induction equation for the purely 
toroidal B field becomes 


— - Vx(vxB) = 

at 


dv r B dv z B 
dr dz 


(31) 


with z the coordinate along the symmetry axis and r the radial 
cylindrical coordinate (Fig. [TJ. The induction equation is iden¬ 


tical to the continuity equation if B is replaced by Ynr (Micono 
et ak j 1999} Bogovalov et aT7||2012)>. Hence, the evolution of B 


along the flow streamline simply follows B oc Ynr, with the pro¬ 
portionality constant set by the initial conditions at the shock. 


3.5. Geometrical factors and orbital motion 

Taking the origin of the coordinate system at the center of the 
massive star, z along the binary axis and r the radial coordi¬ 
nate perpendicular to the binary axis then the flow element is 
at (r cos 0, r sin 9, z) where 9 is the azimuthal angle around the 
binary axis of symmetry (see Fig. [I}. The unit vector e* is 



1.0 1.5 2.0 2.5 3.0 


Figure 2. The RHD flow in the shocked pulsar wind. The star is 
at the origin (0,0) and the pulsar is at (1,0) in units of the or¬ 
bital separation d. A few selected streamlines have been plotted. 
Three different regions have been colored. They correspond to 
the bow shock (blue), the reflected shock (light blue), and the 
back shock (dark blue) regions. 


For each orbital phase, the simulation is scaled with the or¬ 
bital separation d: 


d = 


a( 1 - e 2 ) 


1 - e sin(w - cj p ) 


(34) 


The orbital parameters are those of LS 5039 (McSwain et al. 
2004} |Casares et al.||2005| |2011)l i.e. the semi-major axis a = 
(GMP z lb /4n z ) l,i with P orb = 3.9 days, M = 1.4M 0 + 23 M 0 the 
total mass, e - 0.35 the eccentricity, a> the true anomaly, and 
u) p - 212° the angle at which the binary is at periastron. We 
divide the orbit into 30 phases, each of which requires a (2D) 
calculation of the particle evolution and a (3D) calculation of 
the observed emission. We use 20 cells for 9, 400 for the elec¬ 
tron Lorentz factor y, 20 for the stellar photon energy e*. The 
calculations were parallelized using OpenMP. 


4. Results 


e* = (sin a cos 9, sin a sin 9, cos a), (32) 

with tana = r/z. Taking advantage of the symmetry around 
the binary axis, the speed in the laboratory frame is v = 
(v r cos9, v r sin9, v z ). It is straightforward to see that the boost 
0* (Eq. [25} that applies to the stellar emission seen in the co¬ 
moving frame does not depend on 9: the evolution can be com¬ 
puted using a set of streamlines taken in a plane including the 
binary axis. The unit vector to the observer is 

e 0 bs = (cos i, sin a> sin i, cos co sin i ), (33) 

with i is the inclination of the system and a> the true anomaly of 
the orbit. 

4 When T p » 1 then T as (y — 2 y 2 y i/2 and p ~ {2- y) Y p N p m e c 2 , 
showing that the kinetic energy of the wind is tapped. 


4. 1. The structure of the shocked flow 

The numerical simulation shows the expected double shock 
structure. Numerical diffusivity, induced by our choice of 
Riemann solver, stabilizes the structure despite the presence of 
a strong velocity shear at the interface between the shocked pul¬ 
sar and stellar winds. The diffusivity leads to gradual mixing 
between the winds i.e. numerical spreading of the contact dis¬ 
continuity, quenching the development of the Kelvin-Helmholtz 
(KH) instability (Lamberts et al.|20TT) . Strong KH mixing could 
impact the emission of the region, for instance by reducing the 
Lorentz factor of the flow, and by generating strong turbulence. 
The fluctuation timescale of the interface would be short since 
the flow is relativistic. However, the strong velocity shear is ac¬ 
companied by a strong density contrast between the dense stellar 
wind and the tenuous pulsar wind. The ratio of the KH growth 
timescale to the advection timescale is oc (p 2 /pi) l,2 Av for two 
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fluids of density p\ and P 2 ^ Pi, sheared by a velocity difference 
At' (see appendix in Lamberts et al. 2012 and |Bodo et al.|2004 for 
the growth rate in the relativistic regime). Hence, the KH growth 
is dampened for high density contrasts, such as that expected be¬ 
tween the tenuous highly relativistic pulsar wind and the dense 
stellar wind, making it debatable whether KH-induced mixing is 
dynamically important in gamma-ray binaries on the scales that 
we consider here. Bosch-Ramon et al. (2012 2014[ > find mixing 
occurs mostly on larger scales in their simulations and attribute 
it rather to instabilities triggered by orbital motion. Large, dense 
clumps in the stellar wind could also affect the shock structure 
and variability (Paredes-Fortuny et al. 2015} , although it is un¬ 
clear whether there is enough time for the clumps to grow before 
reaching the termination shock in LS 5039 (located within 1-2 
stellar radii of the star). Our simulation implicitely assumes that 
any mixing is limited and roughly captured — in a time-averaged 
sense — by the numerical diffusivity. 

The basic structure of the shocked pulsar wind is illustrated 
in Fig. [2] where we show only part of the full simulation domain. 
The stellar wind (shocked and unshocked) and the unshocked 
pulsar wind have been edited out of this map as well as sub¬ 
sequent ones since our focus is entirely on the shocked pulsar 
wind. The head of the pulsar wind is shocked in a bow-shaped 
region with asymptotic angles « 40° (termination shock) to 50° 
(contact discontinuity), measured from the z axis. T his is larger 
than the 30° to 45° found by Bogovalov et al. (2008 1 for r\ = 0.1, 
but our simulation domain is smaller and our angles may not yet 
have reached their true asymptotic values. 

Besides the bow-shaped shock, our simulation shows that 
the pulsar is also terminated at the back instead of propagat¬ 
ing freely. The structure is the one expected for Mach reflec¬ 
tion on the binary axis as described in Bogovalov et al. ( 2008} 
in the context of gamma-ray binaries and as observed in e.g. 
simulations of pulsar bow-shock nebulae (Gaensler et al .|2004 ). 
Material flowing in the bow shock region abruptly changes di¬ 
rection when it crosses into the light blue region (black stream¬ 
lines in Fig. [2}. The change is due to a reflection shock that ap¬ 
pears in order to accommodate the back shock region. This re¬ 
flected shock region is separated by a contact discontinuity from 
the back shock region (boundary between medium and dark blue 
regions in Fig.[2|. 

The back shock structure in our simulation is very similar to 
the back shock structure in the 2D and 3D relativistic simula¬ 
tions of |Bosch-Ramon et al.| ( |2012| |2014} , who use a different 
code (PLUTO) but similar values for r\ and F p (rj = 0.1 and 
F p = 2 for their 3D simulation). All their simulations include 
orbital motion and they interpret the presence of this structure as 
an effect of orbital motion. Since this cannot be the case in our 
simulation, we suspect that the confinement depends on a subtle 
combination of 3D + relativistic + pressure (Mach number) ef¬ 
fects ( [Lamberts et al.|2011||2012[|2013} . We defer a resolution 
of this possible issue to future studies. In the present case, as we 
shall see, the back shock and reflected shock only play a minor 
role in shaping the high-energy emission. 

Figure [3] shows maps of the various flow quantities in the 
shocked pulsar wind. The jumps in density single out the re¬ 
flected shock region. The jump in pressure identifies the inter¬ 
face with the bow shock flow as a shock while the matching 
pressures identifies the interface with the back flow as a con¬ 
tact discontinuity. The bow shock flow is re-energized by adia¬ 
batic compression when it crosses the reflected shock. The mag¬ 
netic field distribution is identical regardless of the assumption 
adopted for B at the pulsar termination shock (j |3.4} . The high¬ 
est magnetic field intensities are found at the contact discontinu- 



■ 

9 = 4/3 5/3 


r=l 2 3 4 5 6 7 


Figure 3. Maps of various quantities in the shocked pulsar wind. 
The particle density n , pressure p, and magnetic field are dis¬ 
played on a logarithmic scale ranging down to 10 4 of the max¬ 
imum value (top colorbar). The magnetic field b\ (resp. b 2 ) is 
calculated using Eq. [29] (resp. Eq. [30} at the shock. The adia¬ 
batic index y and Lorentz factor F are displayed on a linear scale 
(bottom colorbars). The map spatial scale is in units of the orbital 
separation d with the pulsar at (1,0) and the star at (0,0). 


ity with the stellar wind, where streamlines from the bow shock 
head pile up. The magnetic field increases with the density in 
the reflected shock region (b oc nr). The last two panels show the 
adiabatic index y and the Lorentz factor. The adiabatic index is 
that of a relativistic gas (y = 4/3) when the shock is perpendic¬ 
ular and decreases towards its non-relativistic value (y = 5/3) 
as material flows in the bow shock region due to adiabatic ex¬ 
pansion. The bow shock flow accelerates back up to a fraction 
of the initial Lorentz factor of the free pulsar wind before be¬ 
ing slowed down again by the reflected shock. The properties in 
the back flow region vary little on the scales examined here: the 
flow speed remains close to v = c/3 with y ~ 4/3, and a slowly 
varying pressure and density. 
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Figure 4. Left: evolution of particle energy as a function of the 
distance / along the streamline (in units of the orbital separa¬ 
tion). Right: evolution of the particle energy distribution along 
the streamline (lighter colors for later times i.e. increasing dis¬ 
tance /). The top panels correspond to the streamline starting at 
o — 48°, the bottom panels to the streamline starting at o = 115° 
(both are identified by an initial dot in Fig. |2j. The distribution 
in the bottom right panel “jumps” when the particles cross the 
reflected shock and are re-energised. The evolution is calculated 
at periastron for our reference simulation. 



Figure 5. Fraction of the mean particle energy remaining after 
radiation losses along each streamline, assuming the conditions 
at periastron ((/> = 0). The streamlines are ordered by angle o 
(Fig# Dashed line shows the same at apastron (</> = 0.5). 


4.2. Particle cooling 


Our reference model for the emission has = 1 (acceleration 
timescale), s = 2 (power-law slope of electron spectrum, no 
Maxwellian), and B calculated using Eq. 29 The parameters C p 
and Cb are set to 1. At the apex of the termination shock at pe¬ 


riastron (d = 0.098 AU), Cp = 1 corresponds to y n 


8 x 1(T 


and Ch = 1 corresponds to B ~ 40 G. Figure [4] illustrates the typ¬ 
ical evolution of particles along two streamlines, at periastron 
<p = 0. We label streamlines by the angle o that their starting 
point makes with the binary axis, with o — 0° corresponding to 
the bow shock head on the binary axis, o = 90° perpendicular to 
the binary axis, 130° < o < 180° to the back shock (see Fig. 2] 
where the black dots identify the two streamlines for which the 
particle evolution is shown). 

For the first streamline (o = 48°), the initial particle distri¬ 
bution ranges from y min = 6.8 x 10 4 to y max = 8.8 x 10 6 (top 
panels of Fig. [4]). The highest energy electrons radiate away half 
of their energy on a scale / ^ 0.003 (in units of the orbital separa¬ 
tion, here d = 0.092 AU). For comparison, this is comparable to 
the spatial resolution at maximum grid refinement in our simula¬ 
tion. Properly resolving the cooling spatial scale within the RHD 
simulation would require an additional 1-2 levels of refinement, 
at large computational cost as mentioned in §3.1 Even the low¬ 
est energy electrons cool on a small scale l £ 0.05 compared to 
the length of the streamline. Synchrotron and inverse Compton 
burnoff at high energy is seen in the evolution of the particle 
distribution (right panel). The late evolution of the distribution 
is set by adiabatic losses, which do not modify the shape of the 
distribution. 

For the second streamline (o = 115°), the initial distribution 
ranges from y min = 1.2x 10 4 to y max = 3.1 x 10 7 . The higher y max 
is due to the lower magnetic field at the shock (Eq.[l9| while the 
lower y m ; n is due to the lower pressure (Eq. 201. Radiative cool¬ 
ing is weakened by the distance and by the higher Lorentz factor 
(decreasing the comoving density of stellar photons). There is 
a moderate evolution of the particle distribution before the par¬ 
ticles are re-energized by passing through the reflection shock. 
Compression at the shock heats the particles and enhances the 
magnetic field. Radiative cooling is much more important in the 
subsequent evolution of the highest energy particles. 

Without radiative losses, the particles lose * 43% of their 
energy adiabatically within the box. The particle energy losses 
are enhanced by radiation. Figure [5 shows the fraction of the 
total energy losses that are due to radiative losses depending on 
streamline. For each streamline, the ratio of the specific energy 
in non-thermal particles e nt to the thermal energy e is compared 
at the beginning and at the end of the streamline. 


6nt (r-1) f r ”’“ 2 dn 

— oc - I ym e c —dy. (35) 

6 P drain ^ 

The value of e nt /e will be the same at the beginning and at the 
end of the streamline (e n t/f)in = (fnt/e)out = Cp if the particle 
losses (or gains) are only due to the adiabatic term i.e. the ra¬ 
tio (e nt /e)i n /(e nt /e)out plotted in Figure [5] will be 1 if there are 
no radiative losses. The figure shows this is not the case in our 
reference simulation: radiative losses dominate the overall en¬ 
ergy losses as the particles propagate along the streamlines. The 
radiative losses are most important for the streamlines that start 
close to the apex (o £ 45°) and in the reflected shock region. 
Integrating e nt N over the whole flow, we find that ~ 70% of the 
power given to non-thermal particles is lost to radiation within 
our simulation box. Adiabatic losses have a minor influence on 
the spectrum and lightcurve: nearly identical results are obtained 
when our baseline calculation is run without taking adiabatic 
losses into account. 

We can only speculate on the feedback that radiative cooling 
could have on the flow dynamics, since our simulation does not 
take it into account. The shock region width is likely to decrease 
as the plasma loses pressure support, raising the density. Since 
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the magnetic field intensity is tied to the density, this may cause 
particles to radiate even faster and at a higher synchrotron fre¬ 
quency than in our computation. Thin shell instabilities may also 
disrupt the interaction region. We leave this for future investiga¬ 
tions. 


4.3. Emission maps 


Emission maps were built for the baseline case with i = 25° 
(Fig. |6}. The maps represent the unabsorbed emissivity inte¬ 
grated over azimuth 6 (Eq. 10 1 


/ 


D 2 ohs j(v/D ohs )rd0. 


This quantity was sampled along each streamline and binned 
to form maps at 1 keV, 1 MeV, 100 GeV, and 1 TeV. Figure [6] 
compares the maps at superior {(f> = 0.08) and inferior conjunc¬ 
tions (0 = 0.77). Animations showing the evolution at all or¬ 
bital phases are available in the online appendix (Figs. A. 1-A.4). 
Fast cooling concentrates emission at the highest frequencies to 
thin layers close to the pulsar termination shock (e.g. compare 
the synchrotron 1 MeV and 1 keV maps). The emission is more 
concentrated at (f> = 0.08 than at 0 = 0.77 because the orbital 
separation is smaller (d ~ 0.10AU compared to d ~ 0.15 AU), 
leading to stronger radiative losses. The simulation box covers 
well the emission zones at 1 TeV and 1 MeV, but misses some 
of the 100 GeV and 1 keV, especially in the back region. Note, 
however, that the map flux scale is logarithmic so the impact on 
the overall lightcurve is negligible. 

Reheating in the reflection shock region is easily seen in the 
maps, especially at 1 keV where a significant fraction of the 
flux may come from this region (and hence escape X-ray ab¬ 
sorption by the stellar wind, see Szostek & Pubus pOTT i. The 
bow shock emission is concentrated towards the head while the 
back shock emission covers a much wider area. The VHE emis¬ 
sion from the back region suffers less from yy absorption and 
actually contributes nearly all the TeV flux when the opacity is 
highest (around <0 = 0.08, see Fig. [8]). 

Figure [7] presents a different way of looking at where 
particles cool. We have plotted the integrated contribution of 
each streamline to the total emission at different frequencies. 
Streamlines that start at o < 90° correspond to the head of 
the bow shock, streamlines with o > 130° correspond to the 
back shock (Fig. [2]). The TeV inverse Compton emission origi¬ 
nates mostly at streamline angles o larger than for the 100 GeV 
emission, and the same applies when comparing the MeV and 
keV synchrotron emission. This is because the streamline initial 
magnetic field decreases with o, allowing for higher initial par¬ 
ticle energies (y rnax ). The back shock contribution clearly dom¬ 
inates the absorbed VHE flux at <0 — 0.08 (dashed lines). At 
<0 = 0.77 the emission is much more concentrated in the stream¬ 
lines with o ~ 100°, which pass through the reflected shock and 
strongly benefit from the relativistic Doppler boost since the flow 
in the back region is then aligned with the observer line-of-sight 
(i = 25°). 


0.0001 0.001 0.01 0.1 

0 = 0.08 (sup. conj.) 0=0.77 (inf. conj. 





1.0 1.5 2.0 2.5 3.0 1.0 1.5 2.0 2.5 3.0 


Figure 6. Emission maps of the shocked pulsar wind at various 
frequencies (top to bottom) and for orbital phases rf> correspond¬ 
ing to the conjunctions (left and right columns). The emission is 
displayed on a logarithmic scale ranging from 1 down to 10 4 , 
with 1 corresponding to the maximum value of the emission 
along the orbit at this frequency (i.e. each map is normalized 
to its maximum value over any pixel and any <0). The VHE maps 
do not take the pair production opacity into account. The case 
shown here corresponds to i = 25° in Fig. [8] The spatial scale is 
in units of the orbital separation d. See Figs. A.1-A.4 in the on¬ 
line appendix for associated movies of the evolution with orbital 
phase. 


4.4. Spectra and lightcurves 

Spectral energy distributions and lightcurves were computed for 
several system inclinations using our baseline model. The re¬ 
sults are displayed in Fig. [8] The spectra and lightcurves are nor¬ 


malised by a coefficient 


7C = 5 x 1(T 10 


p.Skpc'T/ E p 
\ O | \7.6 x 10 35 ergs _1 


erg cm 2 s 1 , 
(36) 
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i = 0° i = 25° i = 50° i = 90° 



eV eV eV eV 



Figure 8. Spectral energy distributions and lightcurves depending on inclination for our reference model. Top panels: thick black 
line is the average spectrum, thick colored lines are the average INFC (dark blue) and SUPC (light blue) spectra (j |4.4| >, thick dashed 
line shows the contribution of the back shock region to the average spectrum, thin grey lines show the spectral evolution with orbital 
phase. Bottom panels: VHE gamma-ray and X-ray lightcurves (thick solid lines). Again, the dashed line shows the back shock 
contribution. The thin grey line in the middle panels is the unabsorbed VHE flux. 


where we have explicited the dependence on the injected power 
in particles E p = 7.6 x 10 35 ergs _1 . E p is related to the pulsar 
spindown power by E = E p ( \ +<x). The results can be scaled with 
E p , or cr, as long as M w and E change in parallel to keep // = 0.1 
(§2). Figure [8] shows the spectral energy distribution sampled at 
various orbital phases to highlight the spectral variability, the av¬ 
erage spectral energy distribution (thick black line) and (in blue) 
the average spectrum corresponding to phases 0.45 < (f> < 0.9 
(INFC) and 0 < 0.45 or <p > 0.9 (SUPC), allowing a comparison 
with the H.E.S.S. spectral analysis in Aharonian et al. (2006). 

The size of the simulation domain limits how far we can fol¬ 
low particle cooling, as the maps of Fig. [6] illustrate. The emis¬ 
sion is thus necessarily incomplete below some energy, which 
we estimate to be Si 100 eV (synchrotron emission component) 
and ^ 1 GeV (inverse Compton component) — based on com¬ 
paring spectra obtained with a reduced domain size. 

The spectra produce broad band X-ray to TeV emission but, 
as could be expected (see |Dubus |2013) ), cannot reproduce the 
peaked GeV emission observed with the Fermi-L AT. This emis¬ 
sion component requires a completely different population of 
electrons, with a narrow distribution in energy. |Zabalza et al.| 
( 2013| ) speculated this could arise from the back shock but we 
find no obvious difference between bow and back shocks. The 
average spectrum from the back region is shown as a dashed line 
in the top panels. Emission from the back region can dominate 
near superior conjunction, when yy absorption is important (see 
dashed lightcurve in the bottom panels), but its contribution to 
the average spectrum remains minor. The spectra of the bow and 
back region are similar ; they would need to have very different 
acceleration parameters to produce significantly different spec¬ 


tra ( |Zabalza et al.|[20l3j |Takata et al.|20l4) . We come back to 
the question of the origin of the HE gamma-ray emission in §5. 

We show lightcurves for the X-ray (1-10 keV) and VHE 
(> 100 GeV) gamma-ray bands, where the spectra and orbital 
modulations are well-known from Suzaku and H.E.S.S. obser¬ 
vations ( |Takahashi et al.||2009[ |Aharonian et al.||2006) l. The 
lightcurves were computed by integrating F(v) over the rele¬ 
vant energy range. When the system is seen face-on (i = 0°), 
the VHE modulation is directly related to the varying stellar 
photon density which increases both inverse Compton emission 
and pair production. The synchrotron emission varies little. The 
synchrotron loss timescale increases at larger orbital separations 
since r sync ~ b -3/2 v l/2 „ d V2 at 

a given frequency. However, the 
actual size of the computational domain increases as d , while the 
particle distribution slope s = 2 ensures equal power per parti¬ 
cle energy, so the emission should vary roughly as d 0 5 , a factor 
1.4 from periastron to apastron, a bit more than what the full 
calculation gives (bottom left panel of Fig. [8j. 

The emission received by the observer changes dramatically 
with the inclination angle of the system. The synchrotron emis¬ 
sion is changed by relativistic boosting as fD 0 b s changes with or¬ 
bital phase. The bow shock region creates a hollow cone of high- 
velocity material, surrounding a filled cone of lower-velocity 
material flowing away from the back shock. Increasing the in¬ 
clination boosts the X-ray emission around inferior conjunction 
0;„r « 0.77, when the shocked flow is oriented towards the ob¬ 
server, and de-boosts the X-ray emission around superior con¬ 
junction </> S up ~ 0.08, when the flow is directed away. At higher 
inclination the observer line-of-sight starts crossing the emission 
cone of the bow shock along its full length. This result in max- 
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o 50 too 150 

o (head: o=0°) 


Figure 7. Contribution to the total flux by each streamline for 
our baseline model with i = 25°. The streamlines are ordered by 
angle o (Fig. [2J. Top panel is for orbital phase rf> = 0.08 (superior 
conjunction), bottom panel for (f = 0.77 (inferior conjunction). 
The flux fractions in each panel correspond to the four frequen¬ 
cies mapped in Fig. [6] Solid lines correspond to the flux unab¬ 
sorbed by pair production, dashed lines is for the absorbed flux 
(affecting the fractions only at 100 GeV and 1 TeV). 


imum boost at rf > m f and 40° ^ i < 50° when the shocked flow 
is going directly in the direction of the observer. For i k, 50° 
the line-of-sight crosses first one edge of the hollow cone then 
the other, resulting in double-peaked emission with a minimum 
at At i = 90°, the line-of-sight crosses the cone first at or¬ 
bital phases 0.54 % <p ^ 0.60 and then at 0.87 $ cf> ^ 0.89, in 
agreement with the position of the two peaks. Although the cone 
is symmetric, the second peak is narrower because of the faster 
orbital motion during the second crossing (nearer to periastron 
passage). 

The VHE emission is influenced by the anisotropy of the 
inverse Compton and pair production cross-sections. Inverse 
Compton emission is enhanced around 0 sup , when stellar pho¬ 
tons are backscattered towards the observer, and diminished 
around cpi„ f, when the stellar photons are forward-scattered 
(Fig# For the same reasons, the pair production opacity is im¬ 
portant around <p sa „ and lower at (/>,„(. The latter can be verified 
by comparing the dark (absorbed) and grey (unabsorbed) lines 
in the middle panels of Fig. [8] However, the effects of Doppler 
boosting dominate the VHE lightcurve. Figure [9] shows the ex¬ 
pected lightcurve at i = 25° but with the relativistic effects turned 
off when computing the emission (the electron populations are 
identical). The X-ray synchrotron lightcurve is nearly constant 
in the absence of Doppler boosting effects. Comparing with the 
full calculation (Fig. [8}, relativistic effects displace the peak 
VHE and X-ray emission towards <t> m i and then lead to a double- 
peaked structure at higher i. Relativistic boosting also hardens 
the VHE spectrum around rf> m f. The intrinsic anisotropic inverse 


i = 25° no boost i = 50° no boost 




0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 1.0 

orbital phase orbital phase 


Figure 9. Same as Fig. [8]for i = 25° and 50° except that relativis¬ 
tic aberrations has not been taken into account when computing 
the emission. 


Compton emission is harder at these orbital phases because scat¬ 
tering is increasingly within the Thomson regime when the stel¬ 
lar photons are closer to being forward-scattered (Dubus et al. 
2008). This effect is amplified by the bulk Doppler boosting. The 


strong spectral evolution with orbital phase at high i can be fol¬ 
lowed in the top panels of Fig# The grey lines show the spec¬ 
tral energy distribution at different phases. The inverse Compton 
spectrum pivots around 100 GeV for i = 50°. Above this energy, 
maximum emission occurs around 0j n f (the INFC spectrum is 
brighter) whereas, below this pivot energy, the gamma-ray emis¬ 
sion peaks around 0 sup (the SUPC spectrum is brighter). The 
inverse Compton lightcurve at different frequencies can thus be¬ 
have in antiphase because of the subtle hardening effects brought 
about by scattering angle and bulk Doppler boost. 


4.5. Exploring parameter space 

We carrried out a limited exploration of the parameter space 
around our reference model. The dependence on inclination i 
has already been shown in Fig. [8] The dependence on the other 
parameters, namely ( p , (/,, £, and s, is shown in Fig. [TO] The 
normalisation of the spectra is the same as for our reference 
model (Eq. 36 1 . We remind that 4 P controls the mean energy of 
the particlesTE q. [20] ), (/, controls the magnetic field intensity at 
the shock (Eq.|29l, controls the maximum particle energy at 
the shock (Eq. 17 1 , and s is the slope of the injected power-law 


distribution of electrons. 

We successively changed the value of each parameter, keep¬ 
ing the others to their reference value. A higher 4b increases syn¬ 
chrotron losses, leading to a pronounced vF v ~ v <2 ~'' ,i2 ~ v° 
spectrum of cooled particles and lowering the inverse Compton 
emission component. Conversely, a lower 4b increases the in¬ 
verse Compton component relative to the synchrotron compo- 
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Figure 10. Dependence of the spectrum on the model parameters. The spectra should be compared to the baseline model with 
— £ — 1,5 = 2 and i = 25° (second column of Fig. [8]). For line labels, see caption to Fig. [8] 


nent and, in our case, leads to a hard synchrotron spectrum 
because the inverse Compton energy losses are in the Klein- 
Nishina regime. A higher £ p increases y mm and thus narrows 
down the energy range of the injected power law. The low- 
energy slope of the synchrotron component for £ p = 10 cor¬ 
responds to the vF v ~ v 4 ' 3 expected for a tail of emission from 
electrons at a high y m j n whereas, in the £ p =0.1 case, the slope 


is the vh\, ~ V ~ v ' slope expected from uncooled elec¬ 
trons emitting in our frequency range with (smaller) Lorentz fac¬ 
tors > y min . Changing the acceleration timescale £ directly im¬ 
pacts the maximum synchrotron frequency, but has no influence 
on the inverse Compton emission because, in our case, y max is 
always high enough for the interaction with stellar photons to 
occur in the inefficient Klein-Nishina regime. Finally, changing 
the slope s of the injected distribution, unsurprisingly, produces 
a harder synchrotron spectrum when the electron distribution is 
harder (smaller s). The ratio of synchrotron to inverse Compton 
emission is also higher because a harder distribution implies 
more very high energy electrons that radiate more efficiently 
synchrotron emission compared to inverse Compton emission in 
the Klein-Nishina regime. 

We do not show how the TeV and X-ray modulations were 
affected by the changes in parameters in Fig. 10 The reason is 
that the modulations did not change much compared to the ref¬ 
erence model. These lightcurves are predominantly shaped by 
the inclination rather than by changes in the other parameters. 
All cases are also comparably radiatively-efficient: the bolomet- 
ric luminosities of the different models vary, in normalised units 
(Eq. [36]), between ~ 13 (£=100) and 35 (£ = 0.01). As in the 
reference case, most of the pulsar power is converted into radia¬ 
tion. 


5. Application to LS 5039 

These results guided us towards a model reproducing the emis¬ 
sion from LS 5039 based on our RHD simulation. This model is 
compared against the observed spectral energy distribution of LS 
5039 in Fig. 11 and the X-ray, MeV, GeV, and TeV lightcurves 
in Fig. 12 The cost of the calculations (several hours per model) 
does not allow an extensive exploration of parameter space. The 


parameter combination given here is only indicative of what 
seems to work i.e. this is not a best-fit model. 


5.1. Model parameters 


The main drivers in deriving this model were (1) reproducing 
the VHE spectral variations ; (2) accounting for the comparable 
X-ray and VHE gamma-ray fluxes ; (3) understanding the origin 
of the HE gamma-ray emission. We start with the latter. 

As the results from the previous section should make clear, 
the HE gamma-ray emission observed with the Fermi-L AT re¬ 
quires an additional component. In principle, a low value of 
£ pushing the synchrotron component to GeV energies might 
account for the Fermi-LAT spectrum with an exponential cut¬ 
off (at the price of supposing a faster-than-Bohm acceleration 
timescale). However, the HE modulation would then be in phase 
with the X-ray modulation, which is ruled out by the observa¬ 
tions. The HE modulation is actually consistent with expecta¬ 


tions for inverse Compton scattering of stellar photons (Abdo 
let al.|2009D . 


We explored the possibility that the HE emission could be 
due to the inverse Compton emission from a narrow Maxwellian 
distribution of electrons, as we had for the case of PSR B1259- 
63 in lDubus & Ceruttili2013l We assumed that a fraction of the 
pulsar wind particles injected at the shock are accelerated to a 
power law, accounting for the broad band X-ray to TeV emis¬ 
sion, while the rest are randomized to this Maxwellian distri¬ 
bution. Adjusting the HE spectrum with the inverse Compton 
emission from the Maxwellian fixes its mean Lorentz factor to 
y t ~ 5000. This also fixes £ /; to ~ 0.017 through Eq. [20] The 
available specific internal energy is a priori identical for the 
Maxwellian and power-law distributions so £ ; , is thus fixed for 
both populations of electrons. The relative contribution of each 
population is adjusted by the fraction of the particle density go¬ 
ing to each (or, equivalently, total energy since e nt is the same). 
The contributions to the average spectrum from each population 
of particles are highlighted in the top panel of Fig. [IT] 

The H.E.S.S. INFC spectrum is hard, best described as a 
power law of photon index of 1.8 combined with an exponen¬ 
tial cutoff at 8.7 TeV. The SUPC spectrum is a steeper power 
law with an index of 2.5 (Fig. [TT|. The two spectra pivot at an 
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energy ~ 200 GeV. Reproducing both the hard INFC spectrum 
and the comparable levels of X-ray and VHE emission turned 
out to be difficult. The models explored in §4 all cut off around 
100 GeV because the high-energy particles responsible for this 
emission are strongly cooled by synchrotron losses. Lowering 
£b decreases the mean magnetic held, hence increases the VHE 
cutoff, but also lowers the X-ray flux relative to the VHE gamma- 
ray flux. There is a trade-off between having enough synchrotron 
emission to account for the X-ray flux, and lowering Q, to enable 
the highest-energy particles to radiate enough VHE photons. As 
discussed in ^4.5| changing £ has little influence on the VHE 
spectrum so we kept £ = 1. Limited exploration showed the best 
agreement was obtained by slightly lowering ^ to 0.5 and by 
taking an injection slope s = 1.5 instead of 2. 


The inclination has an effect on both lightcurves and spectral 
variations. The VHE spectral variations are more pronounced 
with higher i although too high an inclination results in a INFC 
spectrum with a flux that is too low (Fig. |8j. A high inclination 
also results in a pronounced dip of emission at <p ml - = 0.77. The 
H.E.S.S. lightcurve appears double-peaked with a shallow mini¬ 
mum at 0; n f, favoring a model with 25° < i < 50° (Fig. |8j. The 
X-ray modulation is single-peaked at 0; n f, favouring models to¬ 
wards the low end of this range of i. Using i = 35° turned out to 
be a good compromise. 


The “best-adjusted” model shown in Fig. |11|12| has C, p ~ 
0.017, = 0.5, {; = 1, s = 1.5, i = 35° and an injected popu¬ 

lation of particles consisting of a Maxwellian plus a power law. 
An additional parameter is the value of rj = 0.1 that we fixed all 
throughout this study. The hard injection spectral slope s = 1.5 
hints at reconnection rather than Fermi acceleration. Adjusting 
the model to the observations requires that the injected power 
in particles is E p ~ 10 35 erg s -1 . The majority of the particles or 
available power (88%) goes to the Maxwellian at the termination 
shock. Only 12% goes to the particles distributed as a power law. 
However, these power-law particles are more radiatively efficient 
than those in the Maxwellian: about 80% of the power injected 
as a power law ends up radiated away within the simulation do¬ 
main (like the models shown in §4) compared with only 25% of 
the power injected as a Maxwellian. Hence, most of the radiation 
comes from the injected power law but most of the power is in 
the Maxwellian. As a consequence, most of the particles evolve 
adiabatically in this model. 


5.2. Comparison to the observations 


The (rough) adjustment provides a reasonable, albeit imperfect, 
description of the data. The X-ray flux is too low by a factor ~ 2 
in the present model. X-ray emission arising beyond our com¬ 
putational domain might account for this mismatch (§4.4 and 
Fig. [6]). A higher /), would raise the X-ray flux, but lower the 
VHE cutoff in the INFC spectrum to values that would not be 
consistent with the observations. The largest discrepancy is the 
COMPTEL data at MeV energies, recently associated with LS 
5039 ( |Collmar & Zhang||2014| ). Although the spectral slope of 
the model, as well as its evolution with orbital phase (Fig. [I2| , 
are compatible with the observations (the COMPTEL spectrum 
is harder at SUPC phases), the flux is clearly underestimated by 
an order-of-magnitude. A higher would raise the synchrotron 
luminosity but a cooling break in the synchrotron spectrum is 
hard to avoid and this would also be at the expense of the hard 
VHE INFC spectrum. If £), is low then the synchrotron spectrum 
can extend from X-ray to MeV, as in Takahashi et al. ( 2009[ > 
where adiabatic cooling is assumed to dominate over radiative 
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Figure 11. Our best-adjusted model to the observed spectral en¬ 
ergy distribution of LS 5039. Top: the black curve is the av¬ 
erage model spectrum ; the synchrotron (~ 1 eV) and inverse 
Compton (~ 1 GeV) contributions from the Maxwellian popu¬ 
lation of electrons are in dark blue ; the contributions from the 
power-law population are in light blue. Bottom: the thick, dark 
blue curve (resp. light blue curve) represents the INFC (resp. 
SUPC) spectrum, the thin grey curves represent the evolution 
of the SED in orbital phase steps of 1/30. From left to right, 
the data shown are: the orbital maximum and minimum X-ray 
bowties from Suzaku ( |Takahashi et al.||2009 1 , the BATSE data 
points and INTEGRAL bowtie (|Harmon et al.|2004 Hoffmann 


et al. 2009|>, the average COMPTEL data points with bowtie 


(Collmar & Zhang |2014 ), the average Fermi-LAT spectrum with 
data points from 100 MeV to 50 GeV and the best-fit power-law 
with exponential cutoff (Hadasch et al. 2012| >, the H.E.S.S. INFC 


(dark blue bowtie) and SUPC spectra (light blue bowtie) with the 
associated data points from 100 GeV to 30 TeV (Aharonian et al. 
2006). 


cooling, but our model shows the inverse Compton component 
would then be too narrow and luminous (see Fig.fTO}. Our model 
is unlikely to be able to account for the observed MeV flux with¬ 
out additional ingredients (discussed in §6), unless the flux is 
contaminated by diffuse emission or the flux from other sources 
due to the poor angular resolution at these energies. Progress is 
much needed in this difficult observational band. 

The orbital modulations from the model are compared with 
the observations in Figure[l2] The fluxes were calculated in each 
band in the units given by the observations. Because of the mis¬ 
match in X-ray and MeV fluxes, we had to multiply the model 
flux by a factor 2 and 10 (respectively) to obtain a level com¬ 
parable to the observations. The lightcurves are reasonably well 
reproduced. A slightly larger inclination would deepen the VHE 
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minimum at 0 ^ 0.7. The Maxwellian component dominates 
in the HE band, with a modulation in anti-phase with the other 
wavelengths. The simulation output shows that the HE emis¬ 
sion is concentrated at the head of the bow shock and is not 
affected much by relativistic boosting. The HE modulation is 
dominated by the variation with phase of the stellar photon den¬ 
sity and scattering angle, resulting in peak emission at ~ </> sup . 
The Maxwellian component also contributes some flux in the 
10-30 MeV range. A small peak at <p » 0.75 is visible in the 
TeV, MeV, and X-ray lightcurves. This appears to be due to the 
observer line-of-sight grazing the top of the emission cone at this 
orbital phase when i = 35°. The small-scale, stable, flaring struc¬ 


tures observed in the X-ray modulation lightcurve (Kishishita 
et al.|2009 1 might thus be directly related to small structures in 


the shocked flow that are probed when our line-of-sight passes 
through. Emission from the back shock region provides some 
residual TeV flux around <p sup , when emission from the bow 
shock region is strongly absorbed by pair production. This is 
still insufficient to explain the VHE detections. Emission from 
the cascade, initiated when the newly-created e + e pairs are able 
to radiate VHE gamma rays, is very likely to be responsible for 
the residual flux at </> sup . Cerutti et al. (2010) found from their 
study of cascade emission that a good fit required an inclination 
i ~ 40°, consistent with the present model. 

The synchrotron emission from the Maxwellian component 
peaks in the visible band, where it will be difficult to detect 
against the bright V = 11.2 companion O star. This emission is 
modulated, with a lightcurve shape (not shown) similar to the 
MeV modulation. The peak-to-peak amplitude of the V band 
modulation is 0.25 mJy. This translates into a 1.3 mmag mod¬ 


ulation, below the current upper limit of 2 mmag (Sarty et al. 

[20TTT l. 


6. Discussion 

6.1. Influence of the hydrodynamics on the flow emission 

Our motivation for developing this radiative code, based on a rel¬ 
ativistic hydrodynamical simulation, was to obtain a more real¬ 
istic and coherent treatment of the emission geometry, adiabatic 
losses, and Doppler boosting. We discuss these points in turn 
below. 

The simulation shows a complex shock structure, fully con¬ 
taining the pulsar wind, with a reflected shock that re-energizes 
the shocked pulsar wind. In our models, most of the highest- 
energy emission remains concentrated towards the head of the 
bow shock where the electrons cool quickly. Hence, both the 
spectral energy distribution and the modulations are predom¬ 
inantly shaped by the head of the bow shock region. This is 
fortunate as the back shock and reflected shock structures are 
likely to have some dependence on our choice of Mach number 
and could change in the presence of orbital motion, strong mix¬ 
ing, dynamically important radiative losses or magnetic fields, 
etc. It is primarily at lower frequencies, notably in X-rays, that 
these structures contribute significantly, when the particles cool 
more slowly (on larger spatial scales) and/or re-heated to mild 
energies. These conclusions depend on the assumptions that we 
made on what happens at the various shocks, namely that there is 
no difference in particle acceleration between the bow and back 
shock, and that the reflected shock only compresses particles. 
The first assumption is likely to be incorrect at some level be¬ 
cause pulsar winds are not isotropic. Some latitude dependence 
is expected in the pulsar wind due, for instance, to differences 
between the propagation of the high-latitude regions and the 



orbital phase 


id 


Figure 12. Comparison between the LS 5039 model (Fig. 
and the observed lightcurves. The model 10-30 MeV and 1- 
10 keV lightcurves were multiplied by a factor 10 and 2, re¬ 
spectively, to match the observations in Takahashi et al. ( 2009| ; 
Collmar & Zhang ( 2014| >. The grey lightcurve n the 0.1-10 GeV 


and 10-30 MeV panels represents the contribution from the rel¬ 
ativistic Maxwellian component. The VHE and HE gamma-ray 
lightcurves are taken from Aharonian et al. ( 2006| l; Abdo et al. 
(2009), respectively. 


equatorial region (defined by the pulsar rotation axis) where the 
pulsar wind is striped and prone to reconnection. This may be 
manifest in a latitude-dependence of the Lorentz factor of the 
pulsar wind, as seems to be required by models of the Crab pul- 
wind nebula (Bogovalov & Khangoulian||2002 Porth et al. 


|2014|, and/or a dependence of the pulsar wind magnetisation 
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cr with distance (Zabal za et al.|[20l3j |Takata et al.||2014| . A 
latitude-dependent e or B could have important observational 
signatures, even if there is no dramatic change in the structure 
of the flow ( |Vigelius et al.|2007)|Bogovalov et al.|2012| ). On the 
second assumption, having the reflected shock only adiabatically 
compress the particles has some justification because it is usu¬ 
ally found to be difficult to accelerate particles at shock in pair 
plasmas except in special circumstances (very low magnetisa¬ 
tion, shock-induced reconnection of the striped wind). Yet, we 
cannot exclude that particles are re-accelerated to a power law, 
with some influence on the overall emission. We leave the ex¬ 
ploration of these possibilities to future work. 

The models we have explored are radiatively very efficient 
despite the fast flow timescale, the size of the emission re¬ 
gion, and the decreasing magnetic field strength with distance. 
Adiabatic losses play a minor role in the high-energy emission, 
excluding them as the main driver of the X-ray and VHE mod¬ 
ulation in LS 5039 as proposed by Takahashi et al. ( 2009j >. In 
principle, strong radiative losses should be taken into account in 
the dynamics of the flow region. These would result in a nar¬ 
rower and denser shocked flow, experiencing a higher magnetic 
field ( §4.2[ ). However, in order to reproduce the spectral compo¬ 
nent seen with the Fermi- LAT, we have proposed that most of the 
particles in the shocked flow are actually injected in the form of 
a Maxwellian component rather than accelerated to a power-law 
distribution. If the particles in the Maxwellian were only ran¬ 
domized at the shock, their mean Lorentz factor y, corresponds 
to the Lorentz factor l’ ;) of the pulsar wind so Y p « y t = 5000. 
These particles do not cool efficiently in the shocked flow and 
actually dominate the energy budgejj Hence, perhaps counter¬ 
intuitively, the flow remains essentially adiabatic and the as¬ 
sumption of the simulation is verified. 

Doppler boosting has a very strong effect in shaping the 
modulation lightcurves. The geometry is basically a rotating 
cone whose emission is boosted at the phases where its wings 
cross the observer line-of-sight. The result is a double-peaked 
modulation at high inclinations, affecting the synchrotron emis¬ 
sion and the VHE inverse Compton emission. In our models, the 
modulation due to the anisotropy of the inverse Compton cross- 
section is more important than the Doppler modulation only 
at lower energies, when the scattering occurs in the Thomson 
regime. The dependence of the lightcurve shape on inclination 
allows us to constrain i to « 35°. A side effect of Doppler boost¬ 
ing is that it contributes to steepening the VHE emission near 
superior conjunction, erasing the dip around a few 100 GeV in 
the SUPC spectrum due to yy absorption and typically seen in 
previous models (e.g. Dubus et al. 2008[ Khangulyan et al.[2008( 


6.2. Towards more realistic models 

Adjusting to the observations of LS 5039 highlights both the dif¬ 
ficulties and the progress to be expected from the present ap¬ 
proach. Parameters that would be well-suited to the observations 
in one energy band are easily discarded as a result of the com¬ 
prehensive approach taken here, because they fail to reproduce 
the observations in another band. 

We have discussed in §5 the difficulties in reproducing both 
the level of the X-ray flux and the hard VHE spectrum. Part 
of the difficulty may be alleviated if cascade emission is taken 
into account. It would contribute to the VHE flux at all phases, 
not only at superior conjunction, and also to the X-ray emission 
via synchrotron radiation from the pairs ( |Bed narek 2007; Bosch- 
Ramon et al.|2008~ Cerutti et al.|2010]). Combining the present 


Yamaguchi & Takahara 2 012[ >. Our simulation assumed a pulsar 
wind Lorentz factor T /; = 7, not quite high enough to obtain 
ultrarelativistic shock conditions along most of the bow shock. 
A higher F p could enhance the contribution from the wings, al¬ 
though we consider this unlikely given the sharp decrease in flux 
as the shock becomes mostly transverse (Fig. |6|i. Even if a simu¬ 
lation with a higher V p is desirable, we do not expect our results 
to change much. 


5 We have not taken into account the inverse Compton emission from 
the particles in the free pulsar wind, which would contribute to the flux 
in the Fermi -LAT band much like the Maxwellian ffakata et al.12014} . 
Less than 50% of the energy is lost to radiation if T,, = 5000 according 
to Fig. 4 of |Cerutti et al.|{2008[ ). This is an upper limit since their ge¬ 
ometry did not include the back shock, hence the pulsar wind was free 
to propagate to infinity. 


model with a 3D cascade model represents a daunting task. 

The COMPTEL flux level presents a similar challenge to 
models. The modulation, in phase with the X-rays and in anti¬ 
phase with the HE gamma rays, excludes that the COMPTEL 
emission arises from the same electrons responsible for the HE 
gamma-ray emission. It is more natural to attribute it to the ex¬ 
tension of the synchrotron emission ( Takahashi et al.|2009| , yet 
it appears very difficult to achieve this without a high magnetic 
field. As explained above in §6.1, more complex dependencies 
of e or B, motivated by the physics of pulsar winds and their 
termination shock, might be able to simultaneously explain the 
VHE emission and the strong synchrotron emission. 

Our results are based on a single simulation with a given 
//. They should hold qualitatively for other values of r/. The 
strongest impact would certainly be on the value of the inclina¬ 
tion required to adjust the observations since the cone opening 
angle depends directly on //. More subtle effects may appear if rj 
changes along the orbit. This is to be expected at some level in 
LS 5039 because the stellar wind is still in its acceleration phase 
when it encounters the pulsar wind at a distance of one stellar ra¬ 
dius from the surface of the star. Using a beta law for the stellar 
wind velocity, we find that this leads to relative changes of 20% 
in Tj. The opening angle of the cone would be slightly higher at 
periastron than at apastron. Much stronger effects are expected 
in the case of gamma-ray binaries like LS 1+61° 303 and PSR 
B1259-63 where the pulsar wind interacts with a dense equa¬ 
torial outflow from the companion star. Modelling these sys¬ 
tems requires orbit-dependent simulations. Besides the orbital 
phase dependency of rj, such simulations will also be able to ad¬ 
dress the impact of orbital motion on the shape of the interaction 
cone. We expect that the leading arm (the part of the shocked 
pulsar wind that moves into the stellar wind due to orbital mo¬ 
tion) will be compressed and the trailing arm will expand (see 
|Lamberts et al.||2012[ |Bosch-Ramon et aL~||20 14| and references 
therein). The impact on the lightcurves should be limited since 
the emission arises mostly from the innermost, less-affected re¬ 
gions. The back shock in our simulation looks similar to the 
“Coriolis shock” identifie d by |Bosch-Ramon & Barkov ( 201 l| l 
and Bosch-Ramon et al. ( |2012| 2014| > in simulations including 
orbital motion. We speculate that the presence of a back shock 
is not related to the Coriolis force, and note that some of our 
previous simulations and those of Bogovalov et al. | ( 2008j > in¬ 
deed show full confinement without orbital motion, albeit with a 
different back shock geometry (§4.1). Dedicated 2.5D (cylindri¬ 
cal) relativistic simulations would be useful to clearly define the 
conditions for full confinement, the shape of the structure, and 
resolve the issue. 

The parameter Ti, imposes the value of the magnetic field 
at the apex of the bow shock: B ~ 20 G at a distance ~ 3 x 
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10 11 cm from the pulsar and at periastron passage. Taking into 
account the toroidal (beyond the light cylinder radius r\c = 
cP/2n) and dipolar (within light cylinder) nature of the mag¬ 
netic field, the intensity at the pulsar surface is Bo ~ 20 (3 x 


10 11 cm/r L c)(rLc/r„ s ) 3 


1.4 x 10 12 (P/ 0.1 s) 2 G where P is the 


pulsar spin period and for a neutron star radius r ns ~ 10 cm. 
This value of By, for a given P = 0.1 s, is standard for rotation- 
powered gamma-ray pulsars ( Abdo et al.|2013] l. 

The total injected power in particles E p ~ 10 35 erg s' 1 of our 
best model is also standard for pulsars detected in gamma rays 
( Abdo et al.|2013) >. However, combining Eq.[29]for the magnetic 
field with E = E p ( 1 + cr), we find that our model requires 1. 
The pulsar spindown power is equally spread between magnetic 
and kinetic energy. Such a value is much higher than has been 
usually assumed in pulsar wind nebulae, starting with the work 
of Kennel & Coroniti ( 1984) >, though it is not necessarily surpris¬ 
ing since pulsar winds are thought to be launched with very high 
values of cr and to convert the magnetic energy to kinetic energy 


as they propagate (the “cr problem”, see e.g. the reviews by Kirk 
et al.|2009 and Arons 2011). The shock is much closer to the pul 


sar in LS 5039 (« 3 x 10 cm) than in pulsar wind nebulae (0.1 
pc in the Crab nebula) so a higher value of cr is not problematic 
per se. A more worrying issue is that with <r ~ 1 the assumption 
of hydrodynamics breaks down. The higher magnetic field at the 
termination shock means that less of the pulsar wind energy will 
be transferred to the particles. A full RMHD simulation should 
be carried out to investigate whether a substantial magnetisation 
alleviates some of the difficulties we have encountered in repro¬ 
ducing the observations. 

Finally, a pulsar spindown power of E — 2 x 10 35 erg s~* im¬ 
plies a stellar wind mass loss rate M„, = 5x 10 2 M e yr 1 given 
7 / = 0.1 and v w = 2000 km s 1 . This is at the low end of the range 
of estimated M w , even if we take into account that the spindown 
power would have to be increased by a factor at least 3 to account 
for the 1-30 MeV luminosity of w 6 x 10 35 (D/2.5 kpc) 2 erg s _1 
(assuming we could reproduce the peculiar spectral shape). The 
mass loss rate estimated from Her measurements give values in 


the range 2 - 75 X 10 M 0 yr (McSwain et al. 


2004 


Casares 


et al. 2005; Sarty et al.|2011|i. However, wind clumping is known 


to bias this estimator, leading to mass loss rates that can be over¬ 
estimated by a factor ~ 20 for 06.5 stars like LS 5039 (Fullerton 
et al. |2006)>. Indeed, the lack of signatures from X-ray thermal 


emission or absorption in the s tellar wind favours the lower end 


of the range of estimated M w (Szostek & Dubus|2011 


Zabalza 


et al.||20lT|l. Alternatively, tj may be smaller than the value 


we 


assumed. 


7. Conclusions 

We have developped a post-processing radiative code to investi¬ 
gate high-energy non-thermal emission based on relativistic hy- 
drodynamical simulations. Our code includes synchrotron emis¬ 
sion, anisotropic inverse Compton emission, the opacity due to 
pair production at VHE, and takes into account relativistic ef¬ 
fects using the velocity field from the simulation. The particle 
energy distribution is evolved according to the adiabatic losses 
derived from the simulation and radiative losses. Our goal was to 
provide a coherent model of the spectral modulations observed 
from X-ray to VHE gamma rays in the gamma-ray binary LS 
5039. 

(i) The simulation shows a complex shock structure even 
when orbital motion is neglected. The pulsar flow is fully con¬ 
fined by a bow shock and a back shock. The presence of the 


back shock induces a reflected shock in the bow shock region. 
The back shock and reflected shock have a limited impact on the 
overall emission. 

(ii) The VHE emission remains very concentrated towards 
the apex of the bow shock and strongly absorbed by pair pro¬ 
duction at superior conjunction. The back shock contribution 
dominates the VHE flux at superior conjunction but its flux is 
insufficient to explain the H.E.S.S. detection without emission 
from the pair cascade. The back shock thus has a very minor 
influence on the gamma-ray emission from the system. The X- 
ray emission region is much larger, which will help smoothe out 


X-ray absorption signatures from the stellar wind (Szostek & 
|Dubus|2011 1 . 


(iii) Doppler boosting plays the major role in modulating the 
X-ray and VHE emission with orbital phase. Its impact is pre¬ 
dominantly set by the inclination of the system i, with double- 
peaked lightcurves expected at high i. We constrain the inclina¬ 
tion of LS 5039 to i « 35°. 

(iv) There is a tension between the hard VHE spectrum and 
the level of X-ray emission as they require differing intensities 
of the magnetic field. This issue is aggravated by the recent 
COMPTEL detection that, if fully attributed to LS 5039, im¬ 
plies an even stronger synchrotron component (hence higher B) 
and a sharp cutoff between 10 and 100 MeV. These observa¬ 
tions cannot be accommodated in our current model. Possible 
options that may ease the issue include: missing X-ray emission 
from the simulation box, a more intense magnetic field in the re¬ 
gions where radiative cooling is strong (leading to a denser flow 
and a more compressed B), contributions from the pair cascade 
triggered by the absorption of VHE gamma rays, a latitude or 
distance-dependent magnetisation cr or wind Lorentz factor I ' p . 

(v) We attribute the Fermi-LAT emission component to par¬ 
ticles randomized to a Maxwellian distribution at the shock, as 
shown by simulations of particle acceleration at pair-dominated 
shocks. This implies that the Lorentz factor of the wind is 
T p ~ 5000. We find that these particles represent the bulk of 
the power injected in particles, with only 12% going to the elec¬ 
trons accelerated to a power law. Synchrotron emission from 
the Maxwellian population produces a weak (~ 1 mmag) orbital 
modulation of the optical flux. 

(vi) The power-law electrons radiate very efficiently, while 
the “thermal” particles lose energy primarily through adiabatic 
losses. A modest overall injected power of a few 10 33 erg s 1 is 
sufficient to account for the broad band X-ray and TeV emis¬ 
sion. For our choice of ;/ = 0.1 this implies a stellar wind mass 
loss rate of the order of 10 3 M Q yr -1 at the low end of currently 
estimated values. 

(vii) This power, combined with the magnetic field inten¬ 
sity required by our best model, implies a pulsar magnetisation 
cr ^ 1. Such a high value supports the picture that has pulsar 
winds launched with high cr, but fails our assumption of hy¬ 
drodynamics. Relativistic MHD simulations will be required to 
further investigate the issue and, perhaps, resolve some of the 
difficulties encountered in reproducing the observations. 

While gamma-ray binaries may be expected to shed light into 
the processes involved in propagation and termination of pulsar 
winds, we believe that robust conclusions will require the type of 
coherent approach linking dynamical and radiative aspects that 
we have explored here. 
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